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Abstract 

Mathematical models for the description, in a quantitative way, of the 
damages induced on the monuments by the action of specific pollutants 
are often systems of nonlinear, possibly degenerate, parabolic equations. 
Although some the asymptotic properties of the solutions are known, for 
a short window of time, one needs a numerical approximation scheme in 
order to have a quantitative forecast at any time of interest. 

In this paper a fully implicit numerical method is proposed, analyzed 
and numerically tested for parabolic equations of porous media type and 
on a systems of two PDEs that models the sulfation of marble in mon- 
uments. Due to the nonlinear nature of the underlying mathematical 
model, the use of a fixed point scheme is required and every step implies 
the solution of large, locally structured, linear systems. A special effort 
is devoted to the spectral analysis of the relevant matrices and to the 
design of appropriate iterative or multi-iterative solvers, with special at- 
tention to preconditioned Krylov methods and to multigrid procedures. 
Numerical experiments for the validation of the analysis complement this 
contribution. 

1 Introduction 

The problem of monitoring, preserving and, when needed, restoring monuments 
and works of art has become more and more relevant in recent years for the 
conservation of our cultural heritage, after the recognition of the negative ef- 
fects of some pollutants on the monuments. Numerous studies made researchers 
and restorers more and more aware that gaseous pollutants, atmospheric par- 
ticulate matter, and some microorganisms can adversely affect the status of our 
monuments. In order to monitor the cultural heritage and for precisely pro- 
gramming the restoration works, it is of paramount importance to be able to 
accurately assess the status of each monument. Along these lines, quantitative 
methods are emerging and making their way into the practice of preservation 
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and restoration. These have the obvious advantage of allowing fair comparison 
of the state of different monuments, supporting the decision process on what to 
restore, clean, etc and on the relative urgency of each case. 

As an example, consider the "black crusts" that grow on marble surfaces as 
an effect of sulfation of the carbonate stone that is turned into gypsum when 
reacting with SO2 in a moist environment. Since urban concentrations of SO2 
can be nowadays more than 100 times higher than the atmospheric basal values, 
this effect has become very important in the last decades. Sulfation can cause 
permanent damage to the monuments because gypsum crusts can be easily 
eroded by rain or (when located in protected areas) can become unaesthetically 
black including particulate matter from the atmosphere and eventually exfoliate 
|Hay82[ IGKPC891 IBLTROO] . 

A better scheduling of cleaning or deeper restoration can be devised if the 
thickness (and composition) of the crust can be forecast in quantitative way, 
providing a way to compute and thus predict the time evolution of the crust. 
The most common quantitative evaluation of the sulfation phenomena that is 
used in practice consists in assuming that the thickness is directly proportional 
to the length of time of the exposure to the pollutants, with a proportionality 
coefficient obtained by fitting data from a large number of monuments |Lip89| . 
Although this may give an average indication good enough for civil buildings, 
the uniqueness and cultural importance of a work of art calls for a more de- 
tailed analysis, that can take into account the local environment to which the 
monument is exposed. 

A mathematical model of the sulfation of marble based on the chemical 
reactions involved was developed by Natalini and coworkers |ADDN041 IGN05j 
at IAC-CNR (Rome) and tested against experiments, see [GSNF08] , 

It is worthwhile to remark that the mathematical model is able to provide 
new information, which partly contradicts the most common quantitative eval- 
uation methods based on data fitting. In particular (see |GN07| ) the asymptotic 
study of the equations in a one dimensional setting reveals that for large times 
the thickness of the gypsum crust does not grow proportionally to the elapsed 
time as in the Lipfert formula, but proportionally to its square root: the speed 
of growth of the crust is significantly reduced as time goes on. Clearly this 
means that a complete removal of the crust will speed up the damage and calls 
for study of optimal strategies for the periodic partial crust removal. 

However the asymptotic analysis does not give enough information on what 
happens for short times and moreover the study is not yet available for complex 
geometries. For example on a corner stone, SO2 penetrates the marble from 
two sides: how does the crust grow? Does it get rounded? How much? And, 
more importantly, what about the fine particulars of decorations or statues? In 
some cases sulfation caused an almost complete loss of details: can the model 
predict the thickness of the crust there and allow the scheduling of an optimal 
conservation strategy? In order to answer the previous questions, we need a 
numerical method to solve the equations of the model developed by the group 
in Rome (see [ADDN04J). This is a system of two equations, one of which is 
nonlinear of parabolic type. 
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In this paper we generalize and apply novel numerical techniques studied 
in |SDSC| to integrate for long times nonlinear, possibly degenerate, parabolic 
equations like those appearing in the model by [ADDN04 . We wish to point 
out that the techniques developed here have applications that go beyond the 
aforementioned model. For example, in the area of planned conservation, they 
could be adapted to numerically investigate the more complete su lfation mo del 
described in |AFNT07] and the consolidation model presented in [CGN+09] . 

In the literature, degenerate parabolic equations have been discretized mainly 
using explicit or semi-implicit methods, thus avoiding to solve the nonlinear 
equation arising from the elliptic operator. A remarkable class of methods 
arise directly from the so-called non-linear Chernoff formula BP72] for time 
advancement, coupling it with a spatial discretization: for finite differences this 
was started in BBR79] and for finite elements by MNV87]. For example the 
numerical scheme analysed in |ADDN04] for integrating the sulfation model be- 
longs to the class of semi-implicit methods. More recently, another class related 
to the relaxation approximation emerged: such numerical procedures exploit 
high order non-oscillatory methods typical of the discretization of conservation 
laws and their convergence can be proved making use of semigroup arguments 
similar to those relevant for proving the Chernoff formula |CNPS07] . 

In this paper we consider two fully-implicit discretizations in time, thus 
solving a nonlinear system at each time step. In order to fix ideas, consider the 
parabolic equation 

^ = V-(23(«)V«), (1) 

where D(u) is a non-negative differentiable function and denote with Cd{u) the 
elliptic operator on the right hand side. Considering a time discretization such 
that At = t n — i™ -1 and denoting with U the numerical solution, we employ 
the (first order accurate) Implicit Euler scheme 

U(t n ,x) - AtC D (U{t n .x)) = U{t n -\x) (2) 

and the (second order accurate) Crank-Nicholson scheme 

U(t n ,x) - ~£ D (U(t n ,x)) = U{t n ~\x) + ^C D {U{t n -\x)) (3) 

(Note that ^ is also known as the Crandall-Liggett formula, after }CL71j .) 

The computation of U(t n ,x) with ^ or ([3| requires to solve a nonlinear 
equation whose form is determined by the elliptic operator and the nonlinear 
function D{u), but the convergence is guaranteed without restrictions on the 
time step At. Due to the nonlinear nature of the underlying mathematical 
model, the use of a fixed point scheme is required and the choice of the faster 
Newton-like methods implies the solution at every step of large, locally struc- 
tured (in the sense of Tilli, see }Til98| and |SC06| ) linear systems. A special 
effort was devoted in |SDSCj to the spectral analysis of the relevant matrices 
and to the design of appropriate iterative or multi-iterative solvers (see |SC93| ). 
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with special attention to preconditioned Krylov methods and to multigrid pro- 
cedures (see |Gre971 ISaa031 IHac851 ITOS01| and references therein for a general 
treatment of iterative solvers). In this paper we will argue that those meth- 
ods can be extended to the case of systems and perform numerical tests on the 
model of [ADDN04] . 

The paper is organised as follows. In section [2] we recall the results of 
[SDSC on scalar equations, and extend them to the case of a scheme which is 
second order in time. In Section[3]we introduce the sulfation model, the implicit 
numerical schemes and the preconditioners for the linear systems. Both sections 
are complemented by numerical experiments in one and two spatial dimensions. 
Finally in Section [4] we point out some possible developments of this work. 



2 Scalar equations 



In this section, we consider the case of a single equation of the porous media type, 
namely ([!]), where D(u) is a non- negative differentiable function. The parabolic 
equation is of degenerate type whenever D(u) vanishes for some values of u. 
For the convergence analysis of our numerical methods, we will require that 
D(u) is at least continuously differentiable, while the exist ence of solutions is 
guaranteed under the milder assumption of continuity VP 7 . Most applications 
of the porous media equation involve D(u) = u m for some positive m. 

For this particular choice, the followin g self -similar exact solutions have been 



computed by Barenblatt and Pattle (see [V07 



u{t, x) = t 



t a/d 



for t > 0, x e 



(4) 



-i + 



where |x 
singular for t 



= VErX 2 



and a = 



k = 



These solutions are 



d(m-l)+2 ' " 2md ' 

0, but for t > represent important reference cases both for 
the analysis of the solutions of ([T]) and for numerical tests. 

Here, results of [SDSCj on scalar equations are recalled and extended to the 
case of a scheme which is second order in time. They will be used in section 
[3j which deals with a system of two equations, since the linear systems arising 
there include, as a subsystem, those considered in this section. 



2.1 Numerical scheme 

In order to obtain a numerical scheme for approximating the solutions of 
we first discretize the time variable with the Crank-Nicholson form ula gen- 
eralising the simpler case of Implicit Euler that was considered in |SDSC] . We 
point out that searching for high order schemes for equation ([!]) would seem at 
first useless, since the exact solution of the equation are in general continuous 
but not differentiable, so any scheme would converge in theory with order 1. 
However, in practice often the exact solution is piecewise regular, allowing an 
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higher order scheme to converge faster than hrst order and in any case to achieve 
better errors than ^ at a given spatial resolution, even if the theoretical order 
of convergence is not reached. 

We then complete the discretization by considering the points x k = a + kh 
in the spatial domain [a, b], where h = (b — a)/ (JV + 1) and k = 0, . . . , N + l, and 
approximating the one-dimensional Laplacian operator with the usual 3-point 
finite difference formula, i.e. 



d_ 

Ox 



D (u)j+l/2 m\j+l/2 - D ( U )]-l/2 filj-1/2 

— - — r — — + o(l) 



D{u) 0+ i/ 2 {u ]+ i - Uj) - £>(u) j _i /2 (uj - Uj-x) 

(D{u J+1 ) + D(uj))(u j+ i - Uj) - (D(uj) + D(v,j-i))(uj - Uj-Q 

2h 2 



'(1) 
(5) 



In order to write down compactly the equations for the numerical scheme, we 
collect in a vector u™ all the unknown values u" = u n (xj). For example when 
Dirichlct boundary conditions are considered, since Uq and ujv+i are known, u n 
has elements, namely u™, . . . , u 1 ^. 

We denote by tridiag^ , a^, 7&] a square tridiagonal matrix of order N 
with entries @k on the lower diagonal, k = 2, • • • , N, otk on the main diagonal, 
k = 1, • • • , N, and 7fe on the upper diagonal, k = 1, • • • , N — 1. We also denote 
with diag JV [afc] the N x N square diagonal matrix with ak on the k th row. With 
this notation, recalling that ([3]) is a second order approximation, 

u" - u"- 1 = ^^(u-)""^^^-^"- 1 + o{M 2 + h 2 ) (6) 



where 
and 



Ld(u) = tridiagj, [D k _ 1/2 , -D k _ 1/2 - D k+1/2 , D k+1/2 ] (7) 

l) j+i/2 = 5 ' 5 °' ' - - ' 

In two dimensions, on a finite grid composed by the (N + 2) x (N + 2) 
points Xij = a + i/iei + jhe 2 , where a is the lower left corner of the domain, 
h the discretization parameter, e; (I = 1,2) unit vectors along the coordinate 
axis and i,j e N, the matrix Ld( u ) approximating V • (D(u n )Vu n ) is penta- 
diagonal. When considering Dirichlet boundary conditions, adopting the usual 
lexicographic ordering of the unknowns u"j, £r>(u) is a TV 2 x N 2 square matrix 
and nonzero entries can be found only on the main diagonal, on the 1 st and N th 
upper and lower diagonal. 

Finally, we point out that the asymptotic spectral properties of the matrices 
arising from this discretization (Ljj in our case), to a large extent, do not 
depend on the choice of the finite difference formula, but really depend on 
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the Locally Toeplitz structure that in turn arises from operator appearing in 
the PDE f jSC06| ). Thus it should be possible to generalize most of the results 
of the following sections on linear solvers and preconditioning to other spatial 
discretizations techniques, including finite element methods. 



2.2 Newton method 

In order to advance the numerical solution from u n 1 to u n , the nonlinear 
system of equations ^ must be solved at each timestep. We achieve this, by 
iterating with the Newton's method for the function 

F(U) =U " \% Ld ^ U ~ ^ i »<»"- 1 > u "" 1 - uB " 1 ( 8 ) 
The Jacobian of F is 

F / (u)=X N (u) + Y N (u) (9) 

X N (u)=I N -^L D(u) (10) 

Y N (u) = -^T N (u)d\ag%(D' k ) (11) 
Tjv(u) = tridiagf [u fe _! - w fc , u fe _i - 2u k + u k+1 ,u k+1 - u k ] (12) 



where T/v is the same matrix of the first order case (e.g. ( 12 ) in one spatial 
dimension). The only difference in two dimensions is that Tjy is pentadiagonal. 
We observe that the first order scheme based on Implicit Euler gives rise to 

F(u) = u- ^L^u-u"- 1 (13) 

This is the case studied in |SDSCj . Since the Jacobian matrix of F differs from 
F' only for the missing | factors in Xn and Yjv, most of the results proved in 
[SDSCj can be adapted to the present setting. In the following we will thus only 
sketch the proofs. 

Our main result is that the Newton method defined by F, initialised with 
u n,o _ u^-i^ i s convergent under a linear restriction on the timestep. In order 
to prove it, we need the following estimate for the norm of the inverse of J. 

Proposition 2.1. Consider F(u) as defined in ([8]), where u is a sampling 
(at a given time t) of a solution u of ^ with D differentiate and having first 
derivative Lipschitz continuous. If in addition, u is differ entiable with Lipschitz 
continuous first derivative, we have that 

WF'iuy^Kd (14) 

for h sufficiently small and under the additional assumption that At < Cooh for 
some Coo > that does not depend on h. 
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Proof. F' (u) differs from F'{u) for F defined in ( 13 ) only by factors \ appearing 
before any Dk+1/2 term. Since these terms are discarded in the estimates for 
the proof of the analogous result for F, the same proof is valid here. See [SDSC 
for the details. □ 

The following result is a classical tool (see OR70j) for handling the global 
convergence of the Newton procedure. 

Theorem 2.2 (Kantorovich) . Consider the Newton method for approximating 
the zero of a vector function F(u), starting from the initial approximation u^°K 
Under the assumptions that 



F '( u (0)) 

>( U (°)) 



II <0. 



and that 



F(u<°> 

F'{u) -F'(v)|| <7||u-v 

/??77 < -x , 



(15a) 

(15b) 
(15c) 

(16) 



the method is convergent and, in addition, the stationary point of the iterations 
lies in the ball with centre u(°) and radius 

1 - yi - 2p vl 

h 

Theorem 2.3. The Newton method for F(u) defined in ^ for computing u n 
is convergent when initialised with the solution at the previous timestep (i.e. 



u 



u' 



) and for At < Ch, for a positive constant C independent of h. 



Proof. The proof of the same statement for F as given in [SDSC09] can be easily 
adapted to the present case. The technique is to first establish estimates (151 
as follows: 

• P < Ci if At < Cooh by Proposition |2T] 

• r\ < fiCiAt by first applying Proposition |2.1| again and then estimating 



• 7 



< 



F(u"<°) 
D' 



F(u n - 1 ) 



At ||Ld( u ™-i)U 



n-1 1 



0(At) 



00 by direct computation as in (SDSC . 



This implies that condition ( 16 ) can be satisfied when choosing At < Ch for a 
sufficiently small positive constant C that is independent of h. □ 

Remark 2.4. Setting the initial guess with the average between u™ -1 and the 
value given by an Explicit Euler step, like 

„n,0 ,,n-l 1 1 ^ t ,,"-1 
= 27? L c(u-i)U 

does not change the convergence ratio, but in practice one needs less iterations 
to reach a given tollerance. 
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2.3 Iterative methods for the linear system 

Of course the Jacobian matrix F' (u"'( s )) is not explicitly inverted at each New- 
ton step, but instead we compute the (s + l) th Newton iterate by first solving 
the linear system 

j^( u ».W) v W = F(u n 'W) 

for v^ s ) and then setting 

u n,0+i) _ u n,(s) r v («) 

The matrix A N = F'(u n -^) is a square tridiagonal (respectively pentadiag- 
onal) NxN (respectively N 2 x N 2 ) matrix when the domain is one (respectively 
two) dimensional. Its spectral properties are crucial in choosing an appropriate 
solver for the linear system. Given the large dimension of the system, we aim 
at an iterative method with an optimal preconditioner, so that we can compute 
v( s ), on average, in a finite number of iterations. 

Moreover An differs from F' only by the factors \ that were missing in the 
matrices studied in SDSC]. It can thus be shown that An is not symmetric, but 
it is dominated by its symmetric part (An + A~ N )/2, which is in turn essentially 
a weighted laplacian. Rather detailed information on the spectrum of A can be 
gained via the theory of Locally Toeplitz Sequences of [T1198 . In particular, 
when At is chosen proportional to h, the sequence of N x N matrices {hA^} 
obtained for increasing number of grid points is Locally Toeplitz (in the sense of 
|Til98j ) with respect to the pair of functions (D(u(x)),2 — 2cos(s)) defined on 
[a, 6] x [0, 2tt]. Hence (see |SDSC| ) we can expect the GMRES method, which 
is picked due to the asymmetry of An as the main iterative solver, to converge 
in 0(v~N) iterations. 

In order to study a preconditioning strategy, first observe that the sequence 
{Xn} of the symmetric parts of An is also Locally Toeplitz with respect to the 
same generating functions. Next recall that any Locally Toeplitz sequence is also 
Generalized Locally Toeplitz, a class which is an closed under inversion, defined 
in |SC06j . Hence both sequences are also Generalized Locally Toeplitz sequences 
and {X^An} is Generalized Locally Toeplitz with generating function 1 and 
thus the singular values are weakly clustered at the point 1. This is enough to 
guarantee the superlinear convergence of the preconditioned GMRES methods, 
but in |SDSC| we also show that the clustering is strong, proving that Xn is 
an optimal preconditioner for solving a linear system with matrix {Ajv} with 
GMRES, i.e. a given error reduction is reached within a number of iterations 
which is independent on the problem size N. 

Unfortunately there is not a fast direct solver for Xn, so we resort to a 
Multigrid Method (MGM) with a Galerkin approach. The MGM jTOSOlj con- 
sist in constructing a solution of a linear system by composing the action of 
simple iterative schemes (like Jacobi or Gauss-Scidel), that are run on the orig- 
inal system and on smaller systems derived from the first one and called coarse 
grid approximations. 

More precisely, in order to solve a linear system Xu = b in K m , one considers 
a finite sequence of integers mp = m > ra\ > 7712 > • • • > mi > and full-rank 
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matrices P^ +1) € R m -+i xm - 



as 



(called projections) and defines the V-cycle method 

u k+1 = MGM (0,u k ,h) 
with MGM defined recursively as follows: 

Ul (out) :=MGM(i,uf n) ,bi) 



If (j = I) Then Solve(A £ u 



(out) 



Else 1 
2 
3 

4 
5 
6 



b,+i 

u (out) 



(in)> 



P (») r . 

^(i+l) 1 ^ 8 



= P, 



>(i) /( (pW \ 
(i+l) yi (»)v r (i+l)J 

MGM(i+l,0„ j+1) b i+1 ) 

- (-Pi+i) y<+i 



(i) 



Step 1 performs some (z/) iterations of an an iterative method (called pre- 
smoother) for rij-dimensional linear systems that we denoted generically as Si, 
chosen for its error dampening properties, which is often taken in the Jacobi or 
the Gauss-Seidel family. Then, step 2 calculates the residual of the proposed 
solution and steps 3-6 define the recursive coarse grid correction, by projection 
(step 3) of the residual, sub-grid correction (steps 4,5), and interpolation (step 
6). Note that only the smallest system (of level £) is solved exactly, while all the 
others are recursively managed by reduction to low-level system and smoothing. 
For more details and generalisations, see e.g. fTOSOl . 

For differential problems it is natural to construct the coarse grid approxi- 
mations with the Galerkin approach, i.e. by considering a sequence of coarser 
and coarser grids with interpolation operators P(n~ reconstructing values of 
the unknown function on the grid of level I from the smaller set values on the 
coarser grid of level l + l. In this paper it is sufficient to consider linear (bilinear 
in two dimensions) interpolation operators. This is known to give rise to an op- 
timal solver for a weighted laplacian. Moreover in |SDSC| we also observed that 
it is not necessary to bring the MGM to convergence, but applying a single V- 
cycle to the GMRES residual is enough to precondition optimally the GMRES 
method. 



2.4 Numerical tests 

We report here some numerical tests supporting the results of the previous 
sections. For TV ranging from 32 to 1024, we integrate numerically ([!]) with the 
Barenblatt initial data for m — 4 from t — to t — 20/32 in one and two spatial 
dimensions, recording the number of Newton iterations, GMRES iterations and 
the error against the exact solution. 

In figure[l]we plot the quantity ||u 1,s — u 1,s_1 || during the Newton iterations 
s = 1, 2, . . . , 30 for computing the first time step. Different symbols and line 
colours correspond to different mesh sizes on the interval (Ilk) and on a square 
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(a) (b) 

Figure 1: Crank-Nicholson scheme: Newton convergence history in one (a) and 
two (b) spatial dimensions 



* ID: Euler 

+ ID: Crank-Nicholso: 
v 2D: Euler 

* 2D: Crank-Nicholso: 
first order 



Figure 2: Error of the numerical scheme in one and two dimensions, comparing 
with Implicit Euler 

domain (|TJz)) . It is clear that the very good tolerance 10 -6 is reached within a 
reasonable number of iterations: 4 in one spatial dimension and 6 in M 2 . This 
good convergence history is to a large extent not dependent on the chosen value 
for m (see also |SDSUj ). 

Using the Crank-Nicholson scheme, we observe a reduction of the error with 
respect to the first order Euler scheme (see Figure[2|, even if the scheme does not 
converge with the expected order 2. This is due to the presence of singularities in 
the exact solution. The least square fit of the errors obtained, in one dimension, 
with the Crank-Nicholson scheme and represented with plus signa in the figure 
gives that the error decays proportionally to 7V~ 15 . In two dimensions the rate 
of convergence is closer to 1, but the errors are nevertheless lower than those 
obtained with the Implicit Euler scheme. 

Finally, in Figure [3] we plot the number of iterations of the methods for 
linear system. For each value of N, we plot the average number of GMRES 
iterations performed by the algorithm (symbols), while the vertical lines span 
from the minimum to the maximum value recorded during the integration. Panel 
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Crank-Nicholson: GMRES iterations 



O Id: M grid 
X 2d:NxNgrid 



4> 



Crank-Nicholson: PGMRES iterations (1 Vcycle) 

O 'd: N grid, preJacobi 
X 2d: NxN grid, preRBGS 



() () O © <p 



(a) 



(b) 



Figure 3: Crank-Nicholson scheme: GMRES iterations inside the Newton steps, 
(a): no preconditioning (the dashed line is the N 1 / 2 slope), (b): 1 V-cycle as 
preconditioner. 



[3^, compares the number of unpreconditioned GMRES iterations on different 
one- and two-dimensional grids with the VN slope. Note that for the largest 
values of TV, the two-dimensional experiments were not possible due to memory 
limitations on a PC with 8Gb of RAM. We employ MGM as preconditioner, 
choosing one damped Jacobi iteration as a presmoother in the V-cycle in the 
one-dimensional experiments, while for the two-dimensional we chose the more 
efficient Red-Black Gauss-Seidel. Panel [3]d shows that applying one MGM V- 
cycle is optimal (constancy of iterations number for different N) , robust (small 
variance in iterations number) and memory efficient (small number of iterations 
require small amount of storage memory in GMRES). 



3 Marble sulfation 

In this section we consider the model for the marble sulfation problem described 
in [ADDN04 . We describe the main features of the model only briefly, referring 
the reader to the original paper for the details and more comprehensive study 
of the properties of the solutions. [ADDN04] consider the (simplified) chemical 
reaction 

CaC0 3 + S0 2 + ^0 2 + 2H 2 — ► CaS0 4 • 2H 2 + C0 2 . 

to account for the transformation of CaC03 of the marble stone into gypsum 
CaS04 • 2H 2 0, that is triggered in a moist atmosphere by the availability of 
S0 2 at the marble surface and inside the pores of the stone. The two main 
variables of the model are c(t, x) denoting the local concentration of calcium 
carbonate and s(t,x) the local concentration of S0 2 . As the reaction proceeds, 
the calcium carbonate concentration is reduced from the initial value cq, as 
CaC03 is progressively replaced by gypsum. Denoting ip and ip g the porosity 
of the pristine marble and of the gypsum, the model assumes that the porosity 
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(b) 



Figure 4: Sample domains for problem (171 are shown for the ID setting (a) 



and 2D setting (b). The "brick pattern" area represents the marble stone, while 
the dotted area is air. The boundary is drawn with a solid line where Dirichlet 
boundary conditions are applied and with a dotted line where free-flow boundary 
conditions are imposed. 



of the intermediate state is well approximated by linear interpolation 

tp(c) = ip g + {tpo - ip g )— = ac+ 13. 

Co 

The constants a and depend on the porosity of the material involved. The 
model considered in |ADDN04] is described by the following system of PDEs: 



' dip(c)s 



dt 



dc 
dt 



^ v (c)sc + dV ■ M»Vs), 
^(c)sc. 



(17) 



The spatial domain x € 51 in which (171 is set represents a piece of marble 



stone for which at least a portion of the boundary dfl is in contact with the 
polluted atmosphere. In particular dtt is in general split into two parts: one 
represents the outer surface of the marble sample, in contact with the air, and 
the complementary part that separates the portion of the marble object of the 
simulation and the rest of the monument. 

Examples of one and two-dimensional such domains are shown in Figure [4j 
In order to study the formation of a gypsum crust on a flat area of a monument, 
it is sufficient to take a 1-dimensional domain like in the left panel of the figure. 
In order to study a corner-like feature of a work of art, like the edge of a 
monument or a long decoration in relief, one would employ a 2-dimensional 
domain as the one in the right panel of the figure. Obviously more complex 
shapes need 3-dimensional domains representing faithfully the volume occupied 
by the marble. 

Boundary conditions are set by imposing the value of s on the outer boundary 
and by imposing free-flow conditions for s on the inner boundary. In particular, 
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as an example of a one-dimensional setting, we take x € = [0, 1] where x — 
corresponds to the outer boundary of the marble stone, in contact with the 
polluted air, and x = 1 the inner side. The boundary conditions arc illustrated 
in Figure [4j they are of Dirichlet type, imposing s(0,t) on the outer boundary 
and of free-flow type §§(1, t) = on the inner side. 

The parameters m s and m c are fixed by the physical properties of the species 
involved in the reaction and make sure that the mass balance is fulfilled. On the 
other hand a represents the reaction rate and it depends (among other things) on 
the moisture of the air and on the temperature. [ADDN04] describes its central 
role in the analysis of the solutions of the model equations. In particular, if 



u(t,x) is a solution of (17) for a given value of a, then u(t,x) = u{t/a 2 ,x/a) is 
a solution of (17) for a = 1. This observation on one hand plays a fundamental 
role in establishing the long-time asymptotics of the solution and would allow 
to perform the simulations with a = 1 and then rescale the numerical solutions 
appropriately to take the reaction rate into account. However, because of its 
role as fundamental physical parameter of the model, here we prefer to keep a 
explicitly into the equations and perform simulations seeking numerical solution 



of the model in the form (17). Moreover this is beneficial in view of the more 
complete model including a non-constant a has been described in [AF NT07] . 

We consider, as in |ADDN04j a = 0.01, /3 = 0.1, d = 1, m s = 64.06, 
m c = 100.09 and use h = 1/64 or h = 1/128 while varying a from 1 to 10 5 . 

3.1 Discretization 

As in the scalar case, we consider first and second order implicit time dis- 
cretizations, namely the Implicit Euler scheme (Crandall-Liggett formula) and 
the Crank-Nicholson scheme. The general setup for the scheme is the same as 
in the scalar case: we discretize the spatial domain and the elliptic differen- 
tial operator with finite differences, write the time-advancement problem as an 
implicit equation and set up a Newton scheme to solve it. This procedure is 
advantageous if the Newton scheme converges in a reasonable number of itera- 
tions and one can devise an optimal prcconditioner for the linear system that 
has to be solved at each Newton step. 

For the space discretization, we denote x^ — + t;h 6 f2. Approximating 
the elliptic operator along the same lines as in ^ , we consider the second order 
finite difference formula 

d x (ip{c)d x s)\ x , = — 2 (18) 

y(c(^-l/2))(s(Zj) ~ s(Xj-i)) 

This in turn suggests that we employ two staggered grids in the domain f2: 

the grid Xj (j £ N) with the unknowns s" for s{t n ,Xj) and the grid 

(j € N) with the unknowns c" +1 ^ 2 for c(£™, Xj + i/ 2 ). For short, we also denote 

= < ^( c j+i/2)- F° r ease 01 reference, we will denote the two grid also by 
"integer grid" and "half-integer grid" . 
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Wc write explicitly the formulas for the Crank-Nicholson time discretization 
([3]), pointing out that the case of Crandall-Liggett (Implicit Euler) scheme Q 
can be similarly dealt with. Thus we consider the scheme that computes s™ and 
c"+i/2 solving the nonlinear system of equations: 



= F( s )(s",c") = $"s™ + — — C n s n + — dL^s T ' 

2 m r 2 



_-„-l g „-l + &t_ dL , 
2 m r 2 v 



= F( c )(s n ,c") = c" -c™- 1 + ^^-S n c n + 



At a en- 1 „n— 1 
2 m, 



where s™ = [a?, sg, . . .] T , c" = [c™ /2 , c™ /2 , . . .] T and 

Pk+1/2 + Vfc-1/2 

$" = diag fc 



C" = diag fc 
L v n = tridiag fc 

S = dias 



Vk+l/2 C k+l/2 ' Vk-l/2 C k-l/2 



-Pk-l/2i Pk-1/2 + Vfc+1/2) ^^fe+1/2 



(20) 

(21a) 

(21b) 
(21c) 
(21d) 

and 

$ depend on c™ and s", either directly or via the (linear) function ip. It is 
important to note that the matrix L v is defined as of (JtJ) , but it depends 
only on the half of the unknowns of the problem: precisely L v depends on c 
and it multiplies s in formula |20| ). 

The two staggered grids represent a sort of finite difference analogue of the 
approximation with PI (for s{x)) and P0 (for c(x)) conforming finite elements 
considered in [ADDN04 . The results obtained here on preconditioning should 
also be applicable with little modifications in that case too. 
Boundary conditions are imposed considering j 



Note that equations (|20|) are not linear, since also the matrices C, S, L, 

function ip 



1,2, ...,7V in (201 and 



assuming at all time steps a given value for Sq (Dirichlet boundary condition 
at x = 0) and that sn + \ = sat_i (homogeneous Neumann boundary condition 



at x = 1). Thus the expressions of Fj and are modified accordingly 
with respect to those in (20), together with the correspondent elements in the 



Jacobian (22 ) 
u 



The 2N unknowns are collected in a vector u with the ordering 
1 The corresponding sparsity structure of 



91,82, . . . , Sn, Cx/2, . . • , C A r_ 1/ / 2 

the Jacobian matrix is illustrated in Figure [5] For the actual implementation 
it is easier to define the "porous concentration" and set the Dirichlet boundary 
condition as Pj-i/2Sj\- = p So = 1. 

Both the Crank-Nicholson and the Crandall-Liggett formulas give rise to an 
unconditionally stable scheme. Following the results previously established, in 
order to solve the nonlinear problem (20) we set up Newton iterations. To this 



end we need the Jacobian matrix, which is naturally split into four N x N block 
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Figure 5: Sparsity structure of the Jacobian matrix (22) 



as 



J = F' = 



J! 



J! 



u 



The entries (disregarding boundary conditions) are: 



<9F 



j ' k dsi 



[J, 



dF 



(*) 



Cc k+1 / 2 



°*j+l/2 



dF 



0) 

J+l/2 



1/2 



Vj+1/2 + ¥?j-l/2 x At a f j + 1/2 C j + 1/2 + ^-1/2^-1/2 

: ^ fe + 5jk 



dAt 



1 2 h 2 1 "^'-V2^,J'-l + (Pj-1/2 + ^+l/2)4j - </?j + l/ 2 4j + l] , 
<P'j+l/2 S 3 S 3k + ( P' j -l/2 S 3 S 3,k+l 

' 2 

At a (Vj + l/2 C J + l/2 + ¥>j+l/2)5jk + (^-l/2 C J-l/2 + ¥ , i-l/2)*j,fe+l 



2 m 
d At 
2 ~h? 



At a 
2 m 



-<Pj+l/2Cj+l/2- 



At a , 

1 + -2" — (^+i/a c i+i/a + Vj+1/2) 



The sparsity structure of the Jacobian matrix with entries defined in (22 1 is 



(22) 



shown in Figure [5] A more detailed analysis of the matrix will be carried out 
in the next section. 



3.2 Solving the linear system 

At each Newton iterations, we have to solve a linear system with matrix J, 
which is not symmetric and thus we employ GMRES as the main Krylov solver. 
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GMRES iterations - Implicit Euler 



Figure 6: Number (average, min and max) of GMRES iterates per timestep 
with A = 1, for 3 values of N. No preconditioner was used in this test. The 
dashed line indicates the iV 1 / 2 slope. 



We observe that the top-left J| block is given by 



At 



-C 



which is very similar to ^ , except that the identity is replaced by the diagonal 
matrix $ with O(l) entries and the tridiagonal Y term of Q is not present, 
corresponding instead to the third term of the Jj! block. The extra term of Jf, 
involving the diagonal matrix C has entries of order At. Hence we expect J| to 
be spectrally not too different from F' of ^ and thus unpreconditioned GMRES 
iterations count to grow as V~N, which is indeed confirmed in Figure |(3j where we 
plot the average, minimum and maximum number of GMRES iterations needed 
in the case a = 1 and for different values of the number of grid points N. A 
least square fit gives jV - 5217 for the number of iterations. 

In order to devise a preconditioning strategy, given the structure of J", we 
can employ a V-cycle on this block, if we can deal optimally with the rest of the 
matrix. To this end we observe that the lower left block has nonzero entries 
only on two diagonals and these decay as O(At), while the bottom right block 
J£ is the identity matrix plus a diagonal matrix with O(At) entries. 



Theorem 3.1. The upper triangular part of J, 



Jf 







J* 



(23) 



is an optimal preconditioner for J , assuming that the function tp(c) is bounded 
away from 0, or equivalently that /3 > 0. 

Proof. First observe that the diagonal blocks are nonsingular, so that 



p- 







(JcY 
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and the preconditioned system has matrix 



P- 1 J = 



i - wr 1 amr 1 (^) 



WW 



where denotes the null matrix and 1 the identity matrix. 

We now show that all entries of P _1 J are negligible except the diagonal 
ones. In fact J° is diagonal with entries equal to 1 + O(At) and thus its inverse 
has the same property. Since is tridiagonal with entries of 0(aAt), the same 
is true for the lower- left block of P^ 1 J. Gershgorin circles arising from the 
lower half of the matrix are thus centred at 1 in the complex plane and have 
radii decaying as O(At). 

We now turn to consider the upper half of the matrix, where we observe that 
|| (J c s ) (J^r 1 (J s c ) Hoc = O(At) and thus it suffices to show that || (J s s ) _1 is 
bounded to conclude that the Gershgorin circles arising from the upper half 
of the matrix are centred at 1 + O(At) in the complex plane and have radii 
decaying as 0{At). 

To this end, split 

j; = z - w = z(i - z~ x w) 

where Z is the diagonal part, which is 



Z = diag fc (z fc ), 



Zk = <Pk + At — ip k 
m c 



W fc 



where ip k = ^(<fi k +i/2 + fk-1/2) and ip k = l(<p k+1 / 2 Ck+i/2 + ^-1/2^-1/2)- 
Since Z is diagonal, we easily get the estimate 

»^^°ri = s°T^( 1+0 (S))= (£) 

Next observe that Z~ 1 W — A |tridiag fe (^-[(^ fe _ 1/ /2, 0,^+1/2]) and thus 



W = 

II OO 



tridiag fe 



[lfc-l/2,0,</?fc+l/2] 

2p k + £ t Vk + h 2 ^k, 



< max 

k 2ip k 



} hVk + V^ k 



< l-Ch 



for some small positive constant C. Therefore 



1 - z- L w) 



OO ^ 

< V ||^ _1 w|| J < — 



3=0 



Ch 



and 



(Z(l 



w))- 1 = (l-z^wy 1 



0(1) 

□ 
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QMRES iterations with BGS+MCM preeond It loner 



GMRES Iterations with BGS-tVcyele precondiliot 
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N of grid points 
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t t 
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Figure 7: Number (average, min and max) of GMRES iterates per timestep 
with A = 1, for 3 values of N. (a) Block Gauss-Seidel and MGM for the upper 
left block was used as a preconditioner. Blue symbols and lines refer to the 
number of inner MGM iterations, (b) GMRES iterations when performing only 
1 V-cycle of the MGM. 



Remark 3.2. When applying the preconditioner, the block triangular system 
Py = b is solved as 

y c = (J^bc y s = (J*)- 1 ^ - J S cYc) 

where the deponents s and c refer to the upper and, respectively, the lower 
half of the vectors. The previous result shows that the spectrum of P^ 1 J is 
strongly clustered at 1 independently on the discretization parameter h and we 
expect the block-preconditioner P to be optimal. For the whole preconditioner 
to be optimal, however we need an optimal solver for the J| block. However 
Jg is the sum of two diagonal matrices and a tridiagonal matrix which is the 
discretization of a laplacian operator, regularised with the (strictly positive) 
function tp(c(x)), and thus has spectral properties close to those of Xjy studied 



in Section 2.3 As in the scalar MGM (e.g. with 1 damped Jacobi as 

presmoother and a Galerkin approach with linear interpolation) is an optimal 
solver for this block. 

Our preconditioner is thus the Gauss-Seidel preconditioner at block level, 
with MGM on the (s, s) block. The (c, c) block does not need an inner pre- 
conditioner since it is diagonal and can be solved directly. This strategy yields 
an optimal preconditioner, i.e. renders the number of GMRES iterations inde- 
pendent from N, as confirmed by the numerical experiments shown in Figure 
[jj . In the panel [jji we employ MGM driven to convergence as solver for the 
(s, s) block in the preconditioner: GMRES converges in 2-3 iterations, requir- 
ing 10-15 MGM cycles at each iteration. In panel panel [7^, we employ only a 
single MGM V-cycle as inner preconditioner for the (s, s) block: the number 
of GMRES iterations grows slightly (6-8), but this procedure is overall more 
efficient. We observe an impressive series of good features: minimal average 
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computational cost, minimal number of iterations, minimal variance in the lat- 
ter number meaning a strong robustness of the procedure. 

3.3 Simulations and performance of the algorithm 

In Figure [8 we plot some typical curves obtained from the simulations with 
the model (Tj\. Note that for bigger values a, the reaction is faster and a 
boundary layer appears. For a = 10 4 , Figure [9] shows the temporal evolution of 
the two main variables: while S0 2 penetrates deeper and deeper into the stone 
(b), calcium carbonates is substituted by the more porous gypsum in a narrow 
spatial band where the curve c(t, x) presents a boundary layer. Once formed, 
this transition region travels towards the interior of the stone (a). The self- 



similarity of the solutions of (17) under rescaling of the temporal and spatial 
variables mentioned at the beginning of Section [3] implies that the boundary 
layer observed for a = 10 4 will also appear for lower values of a, if the solutions 
were sought for a larger temporal and spatial domain. 



19 



Number of Newton iterations 









O "=i 




i 




# 4=10= 




i 
i 






- * 


i 

i 
i 








! 






- * 








4 




i 





N of grid points 



Figure 10: Number (average, min and max) of Newton iterates per timestep. 



Newton iterations In Figure [10] we plot the average, minimum and maxi- 
mum number of Newton iterations used by the numerical method to solve the 



nonlinear equation (20) at each timestep. 

For a = 1 (black circles) we note that the number of iterations is almost 
constant when the number of grid points is increased. Furthermore, for a given 
number of points employed in the discretization, the number of Newton itera- 
tions increases very moderately even when a is increased by several orders of 
magnitude: e.g. for N = 128 we need an average of 3 Newton iterations per 
timestep for a = 1, 5 for a = 100 (blue stars), and 10 for a = 10000 (red 
crosses). Finally we point out that for higher values of N, the number of New- 
ton iterations decreases slightly since the grid becomes able to resolve better 
the boundary layer. 



3.4 Asymptotics for the front position 



The asymptotic analysis of [GN07j predicts that the front of the travelling wave 
of c(t, x) that separates the gypsum dominated phase from the carbonate dom- 
inated phase and that moves inwards in the marble sample asymptotically be- 
haves as £f r ont ~ v£< 

We performed numerical experiments to test this prediction and to check how 
fast the front approaches this asymptotics. In order to perform the comparison, 
we extracted the information on the front position from the numerical solutions 
c"+i/2 by identifying the gypsum-carbonate front with the point with steepest 
gradient of c(t n , x). 

For a = 10 4 , Fi gurefiT^ shows the position of the front. Note that the step- 
like behaviour of the numerical front that is apparent in some regions of the 
graph is due to the finite spatial resolution of the simulation (h = 1/128). In 
order to check the asymptotics, we plot the position also in double logarithmic 
scale in Figure 11 d, together with the \fi slope (dashed line). We note that both 
simulations agree with the slope of the asymptotics and that for the smaller value 
of a, the solution approaches the asymptotics more slowly. 
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(a) 



(b) 



Figure 11: Gypsum-carbonate front position in a sulfation problem: numerical 
simulations (solid lines) and predicted asymptotics (dashed line) in b 



Number of GMRES and PGMRES(Vcyde) iterations 



5 10 15 20 23 30 33 40 45 

number of time slep 



Figure 12: Number of GMRES iterations per Newton step: for each timestep we 
plot the average, minimum and maximum number of linear iterations. Without 
preconditioner: blue, crosses. With V-cycle preconditioner: black, circles. 



3.5 Sample application in 2D 



In this section we present a numerical simulation of equations (17) in the two 



3.1 



COn- 
TV 

i,3=0 



dimensional setting of figure [4] We consider again two staggered quadrangular 
regular grids in f2 = [0, 1] x [0, 1]. When using N points per direction, we denote 
x £ = €/N and y^ = £/N. We generalize the construction of Section 
sidering two staggered grids: the integer grid is the set of points {(scj, 
carrying the values Sij of the SO2 concentration field s(x, y) and the half-integer 
grid is the set {(x i+ x/ 2 ,yj+i/2)}^j=i carrying the values c l+ i/ 2 j + i/ 2 of the cal- 
cium carbonate concentration field. The discretization of the elliptic operator is 
then generalized in the usual way to the two dimensional setting, and the new 
form of the fixed point problem (20) and its laplacian (22) are derived. The 



numerical scheme now requires, at each timestep, the solution of a system of 
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Carbonate concentration 



marble 




Figure 13: Simulation of marble sulfation in two dimensions. 



2N 2 nonlinear equations. Using again the Newton method, at each iteration we 
need to solve a sparse linear system with a matrix of dimension 2N 2 x 2N 2 . 

The Jacobian matrix has the same block structure described in the one- 
dimensional case (see also Figure [5]). Since it is not symmetric, we use GMRES 
as main Krylov solver with specialised structured preconditioners as in the one 
dimensional case. The main difference in fact is that now the J| block is a 
weighted two-dimensional laplacian plus diagonal corrections that are small (in 
the sense of order of h) . The J£ block remains diagonal with elements equal to 
1 + O(At) and the preconditioner (23 1 can be applied. The off-diagonal blocks 



have more non-zero diagonals, but their elements are still small and Theorem 
|3.1| can be generalized to the two-dimensional setting. 

Here we consider only the best preconditioner of those evaluated in the one 
dimensional setting, namely the upper triangular part of the Jacobian matrix, 
where we perform only 1 V-cycle on the (s, s) block. In Figure 12 we study the 
effectiveness of this preconditioning technique. On a 32 x 32 grid, we observe 
that unpreconditioned GMRES requires an average of 15 to 20 iterations to 
solve the Jacobian linear system in each Newton step, with frequent peaks of 24 
iterations (blue crosses in the figure are the average values, blue lines the min- 
imum to maximum range). Moreover the number of iterations is not constant 
but depends on the timestep. On the contrary the preconditioned method em- 
ploys always an average of 12 PGMRES iterations, with little variability both 
within the time step and across the different times (black circles and lines). 

In Figure [13] we plot the solution obtained for a = 10. We recall the the 
marble is in contact with the polluted air at the bottom and left boundary (cyan 
regions), while at the top and right boundary we apply free flow conditions. 
Both the colour code and the isolines refer to the carbonate concentration in 
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the stone. We observe a clear deformation of the CaC03 field near the corner, 
clearly indicating that SO2, penetrating from both sides, causes an enhanced 
loss of material: if the gypsum crust were to fall off here, the sharp edge would 
be chipped off and the shape of the stone would be permanently changed. This 
simulation, although performed at low resolution (32 x 32 grids) and with a 
moderate value of a, already indicates the relevance of our project of developing 
accurate numerical simulators for realistic geometries of the domain Q in two 
and three dimensions. 

4 Conclusions and future developments 

The novel contribution of this paper relied in the proposal of a fully implicit 
numerical method for dealing with the nonlinear PDE, in its convergence and 
stability analysis, and in the study of the related computational cost. Indeed 
the nonlinear nature of the underlying mathematical model required the appli- 
cation of a fixed point scheme. We have identified the classical Newton method 
in which, at every step, the solution of a large, locally structured, linear system 
has been handled by using specialised iterative or multi-iterative solvers. In 
particular, the spectral analysis of the relevant matrices has been crucial for 
identifying appropriate preconditioned Krylov methods with efficient V-cycle 
preconditioners. Numerical experiments for the validation of our analysis com- 
plement this contribution, which is aimed to provide a non-invasive tool for a 
quantitative forecast of the damage evolution in a given monument. 

In particular we considered the application of the above-mentioned tech- 
niques to the numerical approximation of a mathematical model describing the 
damage of marble monuments by the sulfation process. The use of our resulting 
fast integration algorithms allows to exploit the model and its predictive power 
for the strategy known as planned conservation, that is the novel approach that 
privileges the study and prevention of the damages to delay and optimise the 
actual restoration works. We showed both one dimensional and two dimensional 
numerical simulations, using simple domains. 

For future work, two main directions appear naturally. The first is to finite 
element methods for the space discretization in order to deal with more real- 
istic domains in 2D and 3D, that can model a real architectural item with a 
complicate geometry. In this setting we expect algebraic (linear) MGM or even 
nonlinear multigrid (FAS) to play an important role. The other natural exten- 
sion of the numerical treatment for the sulfation problem involves considering 
the 3-equations model in [AFNT07] and/or a model for a remediation technique, 
like the one in [CGN+09] (this will include consolidation models that is systems 
of the type u t — (D(u)(p(u) x )) x ). As a long term goal, being able to simu- 
late both the damage and the remediation process with validated mathematical 
models and numerical methods would allow to perform numerical experiments 
of restoration works. 
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